Delete Assay 1
- How to remove certain rows for subsequent analysis?
- Maybe just opt to remove them during each analysis or plot using arguments to drop certain rows?
#See change after correction medians align!
par(mar = c(4, 4, .1, .1))
with(D, boxplot(PhiPSII~Assay))
with(D[D$Assay != 1,], boxplot(PhiPSII~Assay))#Dropped Assay 1
# Density plots;transparent
ggplot(D[D$Assay != 1,], aes(x=PhiPSII, fill=Veg.Type)) + geom_density(alpha=.6)#


#————————- # SUMMARY STATS ## Function * New way to enter functions after 2019 packages update:
summary(D)# means, min, max, etc for all variables
## Time Moss Moss.number Moss.sub Microsite
## Min. :0 10 : 1 Min. : 1.00 Length:94 BB-F-37 : 1
## 1st Qu.:0 11 : 1 1st Qu.:22.25 Class :character BB-F-38a: 1
## Median :0 12 : 1 Median :47.50 Mode :character BB-F-39a: 1
## Mean :0 13a : 1 Mean :47.12 BB-F-40a: 1
## 3rd Qu.:0 14a : 1 3rd Qu.:71.75 BB-F-41a: 1
## Max. :0 15a : 1 Max. :96.00 BB-F-42a: 1
## (Other):88 (Other) :88
## Old.Habitat Habitat Index Site Veg.Type
## Length:94 Canopy :48 Min. : 3.00 1:26 Low-scrubland:26
## Class :character Dripline :28 1st Qu.: 8.00 2:34 Mid-shrubland:34
## Mode :character Interspace:18 Median :14.00 3:34 High-woodland:34
## Mean :14.44
## 3rd Qu.:20.00
## Max. :26.00
##
## Topography Shade.Index PDIR Slope Aspect
## S:23 Min. : 6.00 Min. :0.7717 Min. : 0.000 Min. : 0.00
## F:37 1st Qu.:12.00 1st Qu.:0.9142 1st Qu.: 2.000 1st Qu.: 90.75
## N:34 Median :20.00 Median :0.9846 Median : 7.000 Median :200.00
## Mean :17.87 Mean :0.9486 Mean : 7.468 Mean :194.88
## 3rd Qu.:23.00 3rd Qu.:0.9997 3rd Qu.:12.000 3rd Qu.:270.00
## Max. :27.00 Max. :1.0496 Max. :18.000 Max. :357.00
## NA's :4
## VegPercent SolarVeg Camera Pair
## Min. : 0.00 Min. : 0.00 Min. : 1.00 Min. : 1.0
## 1st Qu.: 22.22 1st Qu.: 29.17 1st Qu.: 6.25 1st Qu.: 9.5
## Median : 38.89 Median : 43.75 Median :11.50 Median :19.0
## Mean : 40.66 Mean : 45.26 Mean :11.73 Mean :18.9
## 3rd Qu.: 60.41 3rd Qu.: 62.50 3rd Qu.:16.75 3rd Qu.:28.5
## Max. :100.00 Max. :100.00 Max. :23.00 Max. :37.0
## NA's :72 NA's :55
## Notes.Urgent Assay.Set Assay Clip
## Length:94 Length:94 Min. :1.000 Min. : 1.00
## Class :character Class :character 1st Qu.:1.250 1st Qu.: 6.25
## Mode :character Mode :character Median :2.000 Median :12.00
## Mean :2.468 Mean :12.32
## 3rd Qu.:3.000 3rd Qu.:18.00
## Max. :4.000 Max. :24.00
##
## Fo Fm Fv Fv.Fm
## Min. : 30.00 Min. : 52.0 Min. : 12.00 Min. :0.1300
## 1st Qu.: 79.25 1st Qu.:129.5 1st Qu.: 45.25 1st Qu.:0.3515
## Median : 96.00 Median :183.5 Median : 93.50 Median :0.4770
## Mean : 99.02 Mean :195.0 Mean : 95.97 Mean :0.4546
## 3rd Qu.:118.50 3rd Qu.:248.8 3rd Qu.:132.00 3rd Qu.:0.5697
## Max. :174.00 Max. :419.0 Max. :265.00 Max. :0.7370
##
## Fs Fm.L Fo.L Fv.L
## Min. : 29.00 Min. : 46.0 Min. : 17.00 Min. : 15.00
## 1st Qu.: 76.25 1st Qu.:101.0 1st Qu.: 52.50 1st Qu.: 42.25
## Median : 92.50 Median :137.5 Median : 71.50 Median : 61.00
## Mean : 96.63 Mean :141.8 Mean : 70.89 Mean : 70.86
## 3rd Qu.:118.00 3rd Qu.:175.5 3rd Qu.: 89.25 3rd Qu.: 95.00
## Max. :190.00 Max. :319.0 Max. :149.00 Max. :211.00
##
## Fv.Fm.L PhiPSII qP qNP
## Min. :0.2230 Min. :0.0430 Min. :0.1330 Min. :-0.8570
## 1st Qu.:0.4050 1st Qu.:0.1900 1st Qu.:0.4903 1st Qu.: 0.1765
## Median :0.4905 Median :0.2685 Median :0.5820 Median : 0.4100
## Mean :0.4881 Mean :0.2933 Mean :0.5826 Mean : 0.3651
## 3rd Qu.:0.5780 3rd Qu.:0.3635 3rd Qu.:0.6707 3rd Qu.: 0.5870
## Max. :0.7360 Max. :0.6930 Max. :0.9460 Max. : 0.7970
##
## NPQ Shrubindex Topoindex Siteindex
## Min. :-0.19100 Min. :-1.0000 Min. :-1.000 Min. :-1.00000
## 1st Qu.: 0.08575 1st Qu.: 0.0000 1st Qu.: 0.000 1st Qu.:-1.00000
## Median : 0.32000 Median : 1.0000 Median : 0.000 Median : 0.00000
## Mean : 0.43293 Mean : 0.3191 Mean : 0.117 Mean : 0.08511
## 3rd Qu.: 0.68800 3rd Qu.: 1.0000 3rd Qu.: 1.000 3rd Qu.: 1.00000
## Max. : 2.54200 Max. : 1.0000 Max. : 1.000 Max. : 1.00000
##
## Multibuffer
## Min. :-2.0000
## 1st Qu.: 0.0000
## Median : 1.0000
## Mean : 0.5213
## 3rd Qu.: 1.0000
## Max. : 3.0000
##
#Stats needed for plotting
SE<- function(x) round((sqrt(var(x)/length(x))),3)
Myfunctions = list (
Mean = function(x) round(mean(x),3),
Max= function(x)round(max(x),3),
Min= function(x)round(min(x),3),
Range= function(x) round((max(x)-min(x)),3),
SE= function(x) round((sqrt(var(x)/length(x))),3),
Upper= function(x) round(mean(x)+SE(x),3),
Lower= function(x) round(mean(x)-SE(x),3),
Var = function(x) round(var(x),3),
SD = function(x) round(sd(x),3)
)
Tables
###Stress (Fv/Fm) & Performance (PhiPSII)
#Fv/Fm Tables
#Site table
T.Site.Fvfm<-D %>%
group_by(Site) %>%
summarise_at(vars(Fv.Fm), .funs=Myfunctions)
kable(T.Site.Fvfm, caption = "Fv/Fm Stats by Site", format="pandoc", padding=2)
Fv/Fm Stats by Site
| 1 |
0.440 |
0.602 |
0.183 |
0.419 |
0.021 |
0.461 |
0.419 |
0.011 |
0.105 |
| 2 |
0.429 |
0.737 |
0.130 |
0.607 |
0.027 |
0.456 |
0.402 |
0.026 |
0.160 |
| 3 |
0.491 |
0.686 |
0.144 |
0.542 |
0.028 |
0.519 |
0.463 |
0.026 |
0.163 |
#Topography tables
T.Topo.Fvfm<-D %>%
group_by(Topography) %>%
summarise_at(vars(Fv.Fm), .funs=Myfunctions)
kable(T.Topo.Fvfm, caption = "Fv/Fm Stats by Topography", format="pandoc", padding=2)
Fv/Fm Stats by Topography
| S |
0.432 |
0.656 |
0.20 |
0.456 |
0.026 |
0.458 |
0.406 |
0.016 |
0.125 |
| F |
0.391 |
0.647 |
0.13 |
0.517 |
0.027 |
0.418 |
0.364 |
0.026 |
0.162 |
| N |
0.539 |
0.737 |
0.25 |
0.487 |
0.018 |
0.557 |
0.521 |
0.012 |
0.107 |
#Microhabitat tables
T.Hab.Fvfm<-D %>%
group_by(Habitat) %>%
summarise_at(vars(Fv.Fm), .funs=Myfunctions)
kable(T.Hab.Fvfm, caption = "Fv/Fm Stats by Habitat", format="pandoc", padding=2)
Fv/Fm Stats by Habitat
| Canopy |
0.516 |
0.737 |
0.174 |
0.563 |
0.020 |
0.536 |
0.496 |
0.018 |
0.136 |
| Dripline |
0.421 |
0.647 |
0.130 |
0.517 |
0.028 |
0.449 |
0.393 |
0.021 |
0.146 |
| Interspace |
0.342 |
0.508 |
0.160 |
0.348 |
0.025 |
0.367 |
0.317 |
0.011 |
0.105 |
#Phi-PSII Tables
##Site Table
T.Site.Phi<-D %>%
group_by(Site) %>%
summarise_at(vars(PhiPSII), .funs=Myfunctions)
kable(T.Site.Phi, caption = "PHi-PSII Stats by Site", format="pandoc", padding=2)
PHi-PSII Stats by Site
| 1 |
0.256 |
0.520 |
0.150 |
0.37 |
0.019 |
0.275 |
0.237 |
0.010 |
0.098 |
| 2 |
0.281 |
0.693 |
0.043 |
0.65 |
0.024 |
0.305 |
0.257 |
0.020 |
0.141 |
| 3 |
0.334 |
0.622 |
0.072 |
0.55 |
0.026 |
0.360 |
0.308 |
0.023 |
0.151 |
#Topo table
T.Topo.Phi<-D %>%
group_by(Topography) %>%
summarise_at(vars(PhiPSII), .funs=Myfunctions)
kable(T.Topo.Phi, caption = "PHi-PSII Stats by Topography", format="pandoc", padding=2)
PHi-PSII Stats by Topography
| S |
0.301 |
0.585 |
0.098 |
0.487 |
0.026 |
0.327 |
0.275 |
0.016 |
0.127 |
| F |
0.238 |
0.621 |
0.043 |
0.578 |
0.020 |
0.258 |
0.218 |
0.015 |
0.122 |
| N |
0.348 |
0.693 |
0.072 |
0.621 |
0.024 |
0.372 |
0.324 |
0.019 |
0.139 |
#Microhabitat table
T.Hab.Phi<-D %>%
group_by(Habitat) %>%
summarise_at(vars(PhiPSII), .funs=Myfunctions)
kable(T.Hab.Fvfm, caption = "Phi-PSII Stats by Habitat", format="pandoc", padding=2)
Phi-PSII Stats by Habitat
| Canopy |
0.516 |
0.737 |
0.174 |
0.563 |
0.020 |
0.536 |
0.496 |
0.018 |
0.136 |
| Dripline |
0.421 |
0.647 |
0.130 |
0.517 |
0.028 |
0.449 |
0.393 |
0.021 |
0.146 |
| Interspace |
0.342 |
0.508 |
0.160 |
0.348 |
0.025 |
0.367 |
0.317 |
0.011 |
0.105 |
| #————- |
——— |
– |
|
|
|
|
|
|
|
| # BOXPLOTS & M |
EANS |
|
|
|
|
|
|
|
|
Life zone site boxplots
- Outlier (Phi-PSII) is Moss 34 at Site 2. How could this reading be SO high? Was Fv/Fm high also? Yes: 0.737.
- Question: what dimensions to save for publication? (fig.width=4, fig.height = 1.7)
#Fv/Fm life zone boxplots
PlotA<-ggplot(D, aes(x=(Site), y=Fv.Fm))+
geom_boxplot(outlier.shape=NA, aes(x=factor(Site), fill= (Site)))+
geom_point(aes(fill=(Site)), position = position_dodge2(width=0.3), size=1,
color="black")+
#stat_summary(fun.y=mean, geom="point", shape=20, size=3, color="white")+
geom_point(data = T.Site.Fvfm, aes(x=Site, y=Mean), color="white")+
geom_errorbar(data=T.Site.Fvfm, aes(x=Site, y= Mean, ymin=Lower, ymax=Upper), width=.2, size=0.7, color="white") +
theme_linedraw()+
labs(y='\u2190 Stress (Fv/Fm)', x="Elevation life-zone (site)")+
#scale_x_discrete(labels=c("1" = "Low-scrubland", "2" = "Mid-shrubland", "3"= "High-woodland"))+
scale_x_discrete(labels=c("1" = "Low-site", "2" = "Mid-site", "3"= "High-site"))+
scale_y_continuous(breaks = seq(0, 0.75, 0.1),
expand=c(0,0), limits=c(-0.01,0.75))+
theme(legend.position = "none", text=element_text(size=15), axis.title.x =
element_blank(), panel.grid.minor = element_blank())+
scale_fill_manual(values =c("#E6AB02", "firebrick3", "dodgerblue3"))# panel.grid.major = element_blank(),
#PlotA
#Phi-PSII life zone boxplots
PlotB<-ggplot(D, aes(x=(Site), y=PhiPSII))+
geom_boxplot(outlier.shape=NA, aes(x=factor(Site), fill= (Site)))+
geom_point(aes(fill=(Site)), position = position_dodge2( width=0.3), size=1,
color="black")+
#stat_summary(fun.y=mean, geom="point", shape=20, size=3, color="white")+
geom_point(data = T.Site.Phi, aes(x=Site, y=Mean), color="white")+
geom_errorbar(data=T.Site.Phi, aes(x=Site, y= Mean, ymin=Lower, ymax=Upper), width=.2, size=0.7, color="white") +
theme_linedraw()+
labs(y='Performance (\u03A6PSII)', x="Elevation life-zone (site)")+
scale_x_discrete(labels=c("1" = "Low-site", "2" = "Mid-site", "3"= "High-site"))+
scale_y_continuous(breaks = seq(0, 0.75, 0.1),
expand=c(0,0), limits=c(-0.01,0.75))+
theme(legend.position = "none", text=element_text(size=15), axis.title.x =
element_blank(), panel.grid.minor = element_blank())+
scale_fill_manual(values =c("#E6AB02", "firebrick3", "dodgerblue3"))# panel.grid.major = element_blank(),
#PlotB
#Plot together
PlotsAB <-ggarrange(PlotA, PlotB, labels=c("a", "b"), widths=1:1, ncol=2, nrow=1, font.label= list(size=15)) #Add labels, change text sizes
PlotsAB<- annotate_figure(PlotsAB, top = text_grob("Summer stress assay vs. elevation-life zone sites", size = 15, rot=0, just= "center")) #Add titles
PlotsAB

Topography Boxplots
#Fv/Fm topo boxplots
PlotC<-ggplot(D)+
aes(x=Topography, y = Fv.Fm, fill = Topography) +
geom_boxplot(outlier.shape=NA)+
geom_point(position = position_dodge2( width=0.3), size=1, color="black")+
stat_summary(fun=mean, geom="point", shape=20, size=3, color="white")+
geom_errorbar(data=T.Topo.Fvfm, aes(x=Topography, y= Mean, ymin=Lower, ymax=Upper),
width=.2, size=0.7, color="white") +
scale_y_continuous(name = "\u2190 Stress signal (Fv/Fm)",breaks = seq(0, 0.7, 0.1),
expand=c(0,0), limits=c(-0.01,0.75)) +
scale_x_discrete(name = "Topography (pooled life zones)",
labels=c("S-facing", "Flat", "N-facing")) +
#ggtitle("Summer stress assay (Series 1.1)") +
theme_linedraw() +
theme(text = element_text(size=15), legend.position = "none") +
scale_fill_manual(values = c("#fe9929", "#666666","darkorchid"))+
theme(panel.background = element_rect(fill = "white"), axis.title.x = element_blank(), panel.grid.minor = element_blank())
#PlotC
#Phi topo boxplots
PlotD<-ggplot(D)+
aes(x=Topography, y = PhiPSII, fill = Topography) +
geom_boxplot(outlier.shape=NA)+
geom_point(position = position_dodge2( width=0.3), size=1, color="black")+
stat_summary(fun=mean, geom="point", shape=20, size=3, color="white")+
geom_errorbar(data=T.Topo.Phi, aes(x=Topography, y= Mean, ymin=Lower, ymax=Upper),
width=.2, size=0.7, color="white") +
scale_y_continuous(name = "Performance (\u03A6PSII)",breaks = seq(0, 0.7, 0.1),
expand=c(0,0), limits=c(-0.01,0.75)) +
scale_x_discrete(name = "Topography (pooled life zones)",
labels=c("S-facing", "Flat", "N-facing")) +
#ggtitle("Summer stress assay (Series 1.1)") +
theme_linedraw() +
theme(text = element_text(size=15), legend.position = "none") +
scale_fill_manual(values = c("#fe9929", "#666666","darkorchid"))+
theme(panel.background = element_rect(fill = "white"), axis.title.x = element_blank(), panel.grid.minor = element_blank())
#PlotD
#Plot together
PlotsCD <-ggarrange(PlotC, PlotD, labels=c("c", "d"), widths=1:1, ncol=2, nrow=1, font.label= list(size=15)) #Add labels, change text sizes
PlotsCD<- annotate_figure(PlotsCD, top = text_grob("Summer stress assay vs. Topography zones", size = 15, rot=0, just= "center")) #Add titles
PlotsCD

Microhabitat boxplots
#Fv/fm
PlotE<-ggplot(D)+
aes(x=Habitat, y =Fv.Fm, fill = Habitat) +
geom_boxplot(outlier.shape=NA)+
geom_point(position = position_dodge2( width=0.3), size=1, color="black")+
stat_summary(fun=mean, geom="point", shape=20, size=3, color="white")+
geom_errorbar(data=T.Hab.Fvfm, aes(x=Habitat, y= Mean, ymin=Lower, ymax=Upper),
width=.15, size=0.7, color= "white") +
scale_y_continuous(name = "\u2190 Stress signal (Fv/Fm)",breaks = seq(0, 0.75, 0.1),
expand=c(0,0), limits=c(-0.01,0.75)) +
scale_x_discrete(name = "Microhabitat type") +
# ggtitle("Micro-habitat buffering") +
theme_linedraw() +
theme(text = element_text(size=15), legend.position = "none") +
scale_fill_manual(values = c("darkgreen", "yellowgreen", "tan"))+
theme(panel.background = element_rect(fill = "white"), panel.grid.minor = element_blank(),axis.title.x = element_blank())
#PlotE
#Phi
PlotF<-ggplot(D)+
aes(x=Habitat, y =PhiPSII, fill = Habitat) +
geom_boxplot(outlier.shape=NA)+
geom_point(position = position_dodge2( width=0.3), size=1, color="black")+
stat_summary(fun=mean, geom="point", shape=20, size=3, color="white")+
geom_errorbar(data=T.Hab.Phi, aes(x=Habitat, y= Mean, ymin=Lower, ymax=Upper),
width=.15, size=0.7, color= "white") +
scale_y_continuous(name = "Performance (\u03A6PSII)", breaks = seq(0, 0.75, 0.1),
expand=c(0,0), limits=c(-0.01,0.75)) +
scale_x_discrete(name = "Microhabitat type") +
# ggtitle("Micro-habitat buffering") +
theme_linedraw() +
theme(text = element_text(size=15), legend.position = "none") +
scale_fill_manual(values = c("darkgreen", "yellowgreen", "tan"))+
theme(panel.background = element_rect(fill = "white"), panel.grid.minor = element_blank(),axis.title.x = element_blank())
#PlotF
#Plot together
PlotsEF <-ggarrange(PlotE, PlotF, labels=c("e", "f"), widths=1:1, ncol=2, nrow=1, font.label= list(size=15)) #Add labels, change text sizes
PlotsEF<- annotate_figure(PlotsEF, top = text_grob("Summer stress assay vs. Habitat types", size = 15, rot=0, just= "center")) #Add titles
PlotsEF

Final boxplots: Stress & performance vs 3 scales of buffering
#Plot A-F together
PlotsA_F <-ggarrange(PlotA, PlotB, PlotC, PlotD, PlotE, PlotF, labels=c("a","b","c","d","e", "f"), widths=1:1, ncol=2, nrow=3, font.label= list(size=20)) #Add letters & up the size
PlotsA_F<- annotate_figure(PlotsA_F, top = text_grob("Summer stress & performance vs. 3 scales of habitat buffering", size = 16, rot=0, just= "center")) #Add bottom/top titles
PlotsA_F
#——————— # VARIANCE TESTS
- Goal: Test null of equal variance across categorical levels. Perform batch analysis.
- Test: Fligner-Killeen homogeneity of variance test: stats::fligner.test(). For when groups lack normal distributions.
- Results: Only Fv/Fm ~ Topography is heteroskedastic (for the photosynthesis Anova’s; ignoring vegetation shade metrics for now).
#batch variance tests for life zone sites (Veg.Type)
V.tests<-D %>% select(Fv.Fm, PhiPSII) %>%
map_df(~ tidy(fligner.test(. ~ D$Veg.Type)),.id="Variable")#put results into a tidy df
kable(V.tests, format = "pandoc", caption = "Site Variance: Fligner-Killeen Equal Variance by Site: T0 summer stress (n=94)", digits= c(0, 3, 4))#plot pretty & round output
#Topography
V.tests<-D %>% select(Fv.Fm, PhiPSII) %>%
map_df(~ tidy(fligner.test(. ~ D$Topography)),.id="Variable")#put results into a tidy df
V.tests<-V.tests %>% select(-parameter, -method) #exclude output parameter & method
#plot pretty
kable(V.tests, format = "pandoc", caption = "Topography Variance: Fligner-Killeen Equal Variance by Topozone: T0 summer stress (n=94)", digits= c(0, 3, 4))#round output
#Microhabitat
V.tests<-D %>% select(Fv.Fm, PhiPSII) %>%
map_df(~ tidy(fligner.test(. ~ D$Habitat)),.id="Variable")#put results into a tidy df
V.tests<-V.tests %>% select(-parameter, -method) #exclude output parameter & method
#plot pretty
kable(V.tests, format = "pandoc", caption = "Habitat Variance: Fligner-Killeen Equal Variance by Microhabitat type: T0 summer stress (n=94)", digits= c(0, 3, 4))#round output
Site Variance: Fligner-Killeen Equal Variance by Site: T0 summer stress (n=94)
| Fv.Fm |
5.289 |
0.0711 |
2 |
Fligner-Killeen test of homogeneity of variances |
| PhiPSII |
5.627 |
0.0600 |
2 |
Fligner-Killeen test of homogeneity of variances |
Topography Variance: Fligner-Killeen Equal Variance by Topozone: T0 summer stress (n=94)
| Fv.Fm |
8.398 |
0.0150 |
| PhiPSII |
1.271 |
0.5297 |
Habitat Variance: Fligner-Killeen Equal Variance by Microhabitat type: T0 summer stress (n=94)
| Fv.Fm |
1.498 |
0.4727 |
| PhiPSII |
4.100 |
0.1287 |
#———————— # ONE-WAY ANOVA’s
- Goal: use a family of univariate 1-way ANOVA’s to test for mean differences in photosynthethetic summer stress (Fv/Fm & Phi-PSII) across categorical predictors: elevation-life zone sites, topography zones (aka topozones), & microhabitat types.
- Unbalanced Design Challenge: missing factor levels (South-slope at Site 1 and interspaces at Site 3) prevent a global 3-way ANOVA. Also, all levels have different sample sizes for all three predictors (another type of unbalanced design).
- Heteroskedasticity Challenge: only Fv/Fm by Topozones is heteroskedastic. For this test I will use Welch’s ANOVA if assumption of normality is met. Otherwise, classic ANOVA is used (whether or not normality is perfect).
- P-values: I explicitly avoid adjusting pvalues for this family of tests because this is exploratory analysis using messy ecological data. We want to report all possible patterns of significance!
Classic ANOVAs
#1-way classic ANOVAs - use var.equal = TRUE
W.Aovs<- D %>% select(Fv.Fm, PhiPSII) %>% #select variables
purrr::map_df(~broom::tidy(oneway.test(. ~ Veg.Type, var.equal= TRUE, data= D)),
.id="Variable")#run classic anova
## Multiple parameters; naming those columns num.df, denom.df
## Multiple parameters; naming those columns num.df, denom.df
W.Aovs<-W.Aovs %>% select(-method)
#Extract Cohen's F from batch
CF<-D %>% select(Fv.Fm, PhiPSII) %>%
map_df(., function(x) {
Mod1<- aov(x ~ Veg.Type, data=D)
Out<-sjstats::anova_stats(Mod1)#car pagackage - calculates Cohen's F
Out$cohens.f
})
CF<-as.data.frame(CF) %>% round(., 3)#make df
CF<-t(CF)#transpose
CF<- CF[,1]
W.Aovs<-W.Aovs %>% mutate(Cohens_F = CF)#Add to output
#Extract eta squared (R-squared)
R.sq<-D %>% select(Fv.Fm, PhiPSII) %>%
map_df(., function(x) {
Mod1<- aov(x ~ Veg.Type, data=D)
Out<-sjstats::anova_stats(Mod1)#car pagackage - calculates Cohen's F
Out$etasq
})
R.sq<-as.data.frame(R.sq) %>% round(., 4)#make dataframe
R.sq<-t(R.sq)#transpose
W.Aovs<-W.Aovs %>% mutate(Rsq = R.sq)#Add to output
#Print Pretty Table
kable(W.Aovs, caption = "Classic ANOVAs vs Site - assumes equal variance", digits = c(1, 1, 3, 3, 4, 5, 3), table.attr = "style = \"color: black;\"") %>% kable_styling()
Classic ANOVAs vs Site - assumes equal variance
|
Variable
|
num.df
|
den.df
|
statistic
|
p.value
|
Cohens_F
|
Rsq
|
|
Fv.Fm
|
2
|
91
|
1.614
|
0.2048
|
0.188
|
0.034
|
|
PhiPSII
|
2
|
91
|
2.702
|
0.0724
|
0.244
|
0.056
|
## Multiple Comparisons
Mod<- aov(PhiPSII ~ Veg.Type, data = D, var.equal= TRUE)
TukeyHSD(Mod)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = PhiPSII ~ Veg.Type, data = D, var.equal = TRUE)
##
## $Veg.Type
## diff lwr upr p adj
## Mid-shrubland-Low-scrubland 0.02457240 -0.058961924 0.1081067 0.7635722
## High-woodland-Low-scrubland 0.07795475 -0.005579571 0.1614891 0.0725895
## High-woodland-Mid-shrubland 0.05338235 -0.024383869 0.1311486 0.2360452
Welch’s ANOVAs: Topography
- use stats::oneway.test() with argument var.equal = FALSE, for Welch Correction, an approximate method of Welch (1951), which generalizes the 2-sample Welch test to the case of arbitrarily many samples. Welch’s correction reduces the denominator degrees of freedom in the F-test.
- Question: I read there are no SS’s in Welch’s ANOVA. Why? *pwr::unknown function but the sjstats package must use it to calculate Cohen’s D!
- Multiple comparisons (Games-Howell): South not different from Flat!
#1-way Welch's ANOVAs
W.Aovs<- D %>% select(Fv.Fm) %>% #select variables
purrr::map_df(~broom::tidy(oneway.test(. ~ Topography, var.equal= FALSE, data= D)),
.id="Variable")#run anova
## Multiple parameters; naming those columns num.df, denom.df
W.Aovs<-W.Aovs %>% select(-method)
#Extract R-sq values from summary.lm(Mod)
R.sq<-D %>% select(Fv.Fm) %>%
map_df(., function(x) {
Mod1<- aov(x ~ Topography, data=D)
Out<-summary.lm(Mod1)
Out$r.squared
})
R.sq<-as.data.frame(R.sq) %>% round(., 4)#make dataframe
R.sq<-t(R.sq)#transpose
W.Aovs<-W.Aovs %>% mutate(Rsq = R.sq)#Add to output
#Extract Cohen's F from batch
CF<-D %>% select(Fv.Fm) %>%
map_df(., function(x) {
Mod1<- aov(x ~ Topography, data=D)
Out<-sjstats::anova_stats(Mod1)#car pagackage - calculates Cohen's F
Out$cohens.f
})
CF<-as.data.frame(CF) %>% round(., 3)#make df
CF<-t(CF)#transpose
CF<- CF[,1]
W.Aovs<-W.Aovs %>% mutate(Cohens_F = CF)#Add to output
#Print Pretty Table
kable(W.Aovs, caption = "Welch's ANOVA Fv/Fm vs Toopography - heteroskedastic", digits = c(1, 1, 3, 3, 4, 5, 3), table.attr = "style = \"color: black;\"") %>% kable_styling()
Welch’s ANOVA Fv/Fm vs Toopography - heteroskedastic
|
Variable
|
num.df
|
den.df
|
statistic
|
p.value
|
Rsq
|
Cohens_F
|
|
Fv.Fm
|
2
|
54.745
|
12.098
|
0
|
0.1932
|
0.489
|
#Multiple comparisons - Games-Howell
library(rstatix)
##
## Attaching package: 'rstatix'
## The following objects are masked from 'package:plyr':
##
## desc, mutate
## The following object is masked from 'package:stats':
##
## filter
games_howell_test(data=D, formula= Fv.Fm ~ Topography, conf.level=0.95, detailed=F)
Classic ANOVA: PhiPSII ~ Topography
- use stats::oneway.test() with argument var.equal = TRUE, for classic ANOVA.
- Multiple comparisons result: Phi-PSII differs only between North & Flat!
#1-way Welch's ANOVAs
W.Aovs<- D %>% select(PhiPSII) %>% #select variables
purrr::map_df(~broom::tidy(oneway.test(. ~ Topography, var.equal= TRUE, data= D)),
.id="Variable")#run anova
## Multiple parameters; naming those columns num.df, denom.df
W.Aovs<-W.Aovs %>% select(-method)
#Extract R-sq values from summary.lm(Mod)
R.sq<-D %>% select(PhiPSII) %>%
map_df(., function(x) {
Mod1<- aov(x ~ Topography, data=D)
Out<-summary.lm(Mod1)
Out$r.squared
})
R.sq<-as.data.frame(R.sq) %>% round(., 4)#make dataframe
R.sq<-t(R.sq)#transpose
W.Aovs<-W.Aovs %>% mutate(Rsq = R.sq)#Add to output
#Extract Cohen's F from batch
CF<-D %>% select(PhiPSII) %>%
map_df(., function(x) {
Mod1<- aov(x ~ Topography, data=D)
Out<-sjstats::anova_stats(Mod1)#car pagackage - calculates Cohen's F
Out$cohens.f
})
CF<-as.data.frame(CF) %>% round(., 3)#make df
CF<-t(CF)#transpose
CF<- CF[,1]
W.Aovs<-W.Aovs %>% mutate(Cohens_F = CF)#Add to output
#Print Pretty Table
kable(W.Aovs, caption = "ANOVA: Phi vs Toopography - equal variances", digits = c(1, 1, 3, 3, 4, 5, 3), table.attr = "style = \"color: black;\"") %>% kable_styling()
ANOVA: Phi vs Toopography - equal variances
|
Variable
|
num.df
|
den.df
|
statistic
|
p.value
|
Rsq
|
Cohens_F
|
|
PhiPSII
|
2
|
91
|
6.483
|
0.0023
|
0.1247
|
0.377
|
#Multiple comparisons
Mod<- aov(PhiPSII ~ Topography, data = D, var.equal= TRUE)
TukeyHSD(Mod)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = PhiPSII ~ Topography, data = D, var.equal = TRUE)
##
## $Topography
## diff lwr upr p adj
## F-S -0.06330905 -0.14529288 0.01867479 0.1624954
## N-S 0.04706138 -0.03629741 0.13042018 0.3741169
## N-F 0.11037043 0.03701940 0.18372146 0.0015646
Habitat ANOVAs
- Equal variances, so classic anovas:
- Multiple comparisons (Fv/Fm): Canopy differs from dripline & interspace (interspace = dripline)!
- Multiple comparisons (Phi): All three habitat types differ!!
#1-way classic ANOVAs - use var.equal = TRUE
W.Aovs<- D %>% select(Fv.Fm, PhiPSII) %>% #select variables
purrr::map_df(~broom::tidy(oneway.test(. ~ Habitat, var.equal= TRUE, data= D)),
.id="Variable")#run classic anova
## Multiple parameters; naming those columns num.df, denom.df
## Multiple parameters; naming those columns num.df, denom.df
W.Aovs<-W.Aovs %>% select(-method)
#Extract Cohen's F from batch
CF<-D %>% select(Fv.Fm, PhiPSII) %>%
map_df(., function(x) {
Mod1<- aov(x ~ Habitat, data=D)
Out<-sjstats::anova_stats(Mod1)#car pagackage - calculates Cohen's F
Out$cohens.f
})
CF<-as.data.frame(CF) %>% round(., 3)#make df
CF<-t(CF)#transpose
CF<- CF[,1]
W.Aovs<-W.Aovs %>% mutate(Cohens_F = CF)#Add to output
#Extract eta squared (R-squared)
R.sq<-D %>% select(Fv.Fm, PhiPSII) %>%
map_df(., function(x) {
Mod1<- aov(x ~ Habitat, data=D)
Out<-sjstats::anova_stats(Mod1)#car pagackage - calculates Cohen's F
Out$etasq
})
R.sq<-as.data.frame(R.sq) %>% round(., 4)#make dataframe
R.sq<-t(R.sq)#transpose
W.Aovs<-W.Aovs %>% mutate(Rsq = R.sq)#Add to output
#Print Pretty Table
kable(W.Aovs, caption = "Classic ANOVAs vs Microhabitat type - equal variance", digits = c(3, 3, 3, 3, 4, 5, 4), table.attr = "style = \"color: black;\"") %>% kable_styling()
Classic ANOVAs vs Microhabitat type - equal variance
|
Variable
|
num.df
|
den.df
|
statistic
|
p.value
|
Cohens_F
|
Rsq
|
|
Fv.Fm
|
2
|
91
|
12.433
|
0
|
0.523
|
0.215
|
|
PhiPSII
|
2
|
91
|
14.078
|
0
|
0.556
|
0.236
|
#Multiple comparisons (Fv/Fm)
Mod<- aov(Fv.Fm ~ Habitat, data = D, var.equal= TRUE)
TukeyHSD(Mod)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Fv.Fm ~ Habitat, data = D, var.equal = TRUE)
##
## $Habitat
## diff lwr upr p adj
## Dripline-Canopy -0.09583631 -0.1715138 -0.02015878 0.0091863
## Interspace-Canopy -0.17431250 -0.2622704 -0.08635459 0.0000250
## Interspace-Dripline -0.07847619 -0.1746205 0.01766808 0.1321162
#Multiple comparisons (Phi)
Mod<- aov(PhiPSII ~ Habitat, data = D, var.equal= TRUE)
TukeyHSD(Mod)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = PhiPSII ~ Habitat, data = D, var.equal = TRUE)
##
## $Habitat
## diff lwr upr p adj
## Dripline-Canopy -0.11274107 -0.1813238 -0.04415837 0.0005046
## Interspace-Canopy -0.15511806 -0.2348298 -0.07540627 0.0000349
## Interspace-Dripline -0.04237698 -0.1295077 0.04475369 0.4807655